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FIELD OF THE INVENTION 

The present invention relates to a method of constraining stochastic models 
representing heterogeneous undergroimd zones such as oil reservoirs to data referred to 
as dynamic because they vary with the fluid displacements. These data are, for example, 
5 production data or pressure data obtained from well tests. 

BACKGROUND OF THE INVENTION 

The use of stochastic models of Gaussian type for representing the heterogeneity of 
underground structures is for example described by : 

- Joumel, A.G. and Huijbregts, Ch. J. : "Mining geostatistics". Academic Press, 1978, 
10 or 

- Chiles, J. P. and Delfiner, P. : "Geostatistics - Modeling Spatial Uncertainty", Wiley- 
Interscience Publishers, John Wiley & Sons, 1999. 

A numerical reservoir model can be formed from a set of grid cells to which the 
values of a realization of a stochastic model of Gaussian or related type are assigned. 
1 5 These values can be assimilated to porosities or permeabilities. 

Matching the numerical reservoir model with the dynamic data measured in the 
field can be done in form of an optimization problem. A previously defined objective or 
cost frmction quantifies the difference between the dynamic data measured for the real 
medium and the corresponding responses of the numerical reservoir model. These 
20 responses are calculated by means of a numerical flow simulator. The goal of the 
optimization problem is to modify the reservoir model or rather the associated 
realization to minimize the objective function. This process is iterative : each iteration 



implies direct simulation of the flows. A good optimization method should allow to : a) 
modify realiaations discretized on a very large number of grid cells; b) carry out the 
modifications while respecting the stochastic model, i.e. the modified realization has to 
be coherent with the stochastic model; and c) limit the number of direct flow 
5 simulations because they require a considerable calculating time. 

Simulated annealing can be mentioned as an example of known optimization 
techniques. This approach is for example described by : 

- Gupta, A. D. et al. : "Detailed characterization of fractured limestone formation using 
stochastic inverse approaches", SPE Ninth Symposium, 1994. 

10 This technique is based on realization values exchanges between grid cells. Upon 

each exchange, the objective function has to be calculated and therefore a direct flow 
simulation has to be carried out. This process requires an excessive number of 
iterations. Furthermore, in order to preserve the agreement between the realization and 
the stochastic model, an additional term concerning the variogram is introduced in the 

1 5 objective fonction, which makes optimization more delicate. 

Other optimization techniques, more commonly applied, are based on gradients 
calcvdation. Several approaches based on gradients are presented by : 

- Tarantola, A. : "Inverse problem theory - Methods for data fitting and model 
parameter estimation", Elsevier Science Publishers, 1987. 

20 They require calculating the gradients of the objective fvmction with respect to the 

parameters of the problem which are the values of the realization at each grid cell. The 
realizations are then modified as a function of these gradients so that the objective 
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function decreases. The problem related to conditioning of a reservoir model to 
production data is not linear : the minimization techniques using gradient calculation are 
used iteratively. After each modification of the realization, a direct flow simulation is 
carried out and the gradients are recalculated. Applied suddenly, the gradient methods 
5 lead to calibration of the dynamic data but they destroy the coherence between the 
stochastic model and tiie realization. Besides, they do not allow to consider a very large 
number of parameters. Li order to overcome these limits, geostatistical parameterization 
techniques can be integrated thereto. The pilot point method can be mentioned at this 
stage, which is described by : 

10 -de Marsily, G. et al. : "hiterpretation of interference tests in a well field using 

geostatistical techniques to fit the permeability distribution in a reservoir model" in 
Verly, G. et al. (ed.), Geostatistics for natural resources characterization, Part 2, D. 
Reidel Pub. Co, 1984. 

This method consists in selecting in the realization a certain number of points 
15 referred to as pilot points, in calculating the derivatives of the objective function with 
respect to the values at these points, in modifying the values of these points accordingly 
and in propagating the disturbance thus defined by means of a kriging technique. The 
pilot point method can induce deviant value variations of the pilot points. 

Another geostatistical parameterization technique, which allows the above- 
20 mentioned difficulty to be overcome, is the method of gradual deformation of a 
stochastic model of a heterogeneous meditim such as an underground zone. It is 
described and used by Hu, L.-Y. et al., in patents FR-2,780,798 and FR-2,795,841 filed 
by the applicant. 



The gradual deformation method allows to gradually modify a realization of a 
stochastic model of Gaussian or related type while respecting this model. The deformed 
realization still is a realization of the stochastic model. When the gradual deformation 
method is introduced in an optimization process, the procedure is as follows. The initial 
5 realization is combined with a fixed number of independent realizations related to the 
same stochastic model. These realizations are called complementary realizations. 
Combination is controlled by as many deformation parameters as there are 
complementary realizations. It produces a new realization. The derivatives of the 
objective function with respect to the deformation parameters are then calculated. The 

10 latter are modified so £is to take into accoimt the information from the derivatives. A 
first optimization in relation to the deformation parameters provides a realization 
verifying the stochastic model and reducing the objective function. In general, this 
optimization process has to be repeated several times with different complementary 
realizations so as to sufficiently reduce the objective function, vsiiich may require, in 

1 5 some cases, a prohibitive number of direct flow simulations. 

The gradual deformation method allows, in some cases, to modify a realization 
locally. This possibility is justified when the gradual deformation is combined with the 
FFTMA geostatistical generator described by : 

- Le Ravalec, M. et al. : The FFT moving average (FFT-MA) generator : an efficient 
20 numerical method for generating and conditioning Gaussian simulations. Math. Geol., 
32(6), 2000. 

This generator produces realizations for a stochastic model of Gaussian type 
specified beforehand by convoluting a Gaussian white noise with an operator depending 
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on the covariance function. A local deformation can be carried out by applying the 
gradual deformation method to the Gaussian white noise imderlying the realization. 

The gradient techniques developed to date for calibration in relation to dynamic 
data are based on a direct link between the variations to which the realization 
5 representing the reservoir is subjected and the variation o^ the objective function. The 
optimization process involves modifying first the realization, then starting a flow 
simvdation to apprehend the resulting variations for calibration. 

A different approach valid in cases where the flows are modelled by streamlines, is 
proposed by : 

10 - Wang, Y. and Kovscek, A. R. : A streamline approach for history-matching 
production data, SPE/DOE TOR, 2000. 

Inverting the procedure of conventional approaches, these authors propose a 
methodology based on the successive application of two st£^es. The first stage focuses 
on calibration : it allows to evaluate the modification to be applied to the efifective 
15 permeabilities of the streamlines in order to improve calibration. The second stage 
relates to the transfer of the effective permeability variation of the streamlines to the 
realization. This process is continued until a satisfactory calibration is obtained. The 
main limit of this approach es that it does not preserve the coherence of the realization in 
relation to the stochastic model. 

20 SUMMARY OF THE INVENTION 

The method according to the invention allows to work out numerical models, 
representing heterogeneous underground media such as oil reservoirs or aquifers, in 
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accordance with a set of dynamic data measured in production, injection or observation 
wells, and with a stochastic model characterizing the spatial distribution of the 
heterogeneities. 

It essentially comprises minimization of the objective function from an iterative 
5 process subdivided in two stages and by integration of the gradual deformation method 
in the second stage. The first stage identifies the disturbance to be applied to the 
effective permeabilities of predetermined zones so as to reduce the objective function. 
The second stage focuses on the minimization of the difference between the desired 
effective permeabilities and the corresponding effective permeabilities calculated for the 
10 realization considered, by modifying this realization by means of the gradual 
deformation method. This minimization stage is an intermediate optimization which 
requires no new flow simulation and allows to deform the realization while respecting 
the Gaussian stochastic model. If the objective function is not satisfactory at this stage, 
the process is resumed at the first stage. 

15 The method according to the invention is a decoupled gradual optimization method. 

The term gradual applies insofar as deformation of the nimierical model is carried out 
according to the fundamental principles of the gradual deformation method. The term 
decoupled reminds that the optinaization is based on two stages, the first one defining 
the correction to be brought to the effective permeabilities, the second propagating it 

20 independently to the numerical model. These two stages imply that the flow simulations 
are carried out only during the first stage in order to improve calibration of the effective 
permeabilities. On no account does the second stage require a flow simulation to 
estimate the way the numerical reservoir model should be deformed. 



In other words, the method according to the invention allows to form iteratively a 
model representative of the permeability field of a heterogeneous medium, discretized 
by a grid, such as an underground zone, constrained by a priori geologic data and 
dynamic data collected in the medium by means of measurements and observations 
5 obtained beforehand. It includes : 

- a first stage comprising : 

a) generating an initial model of the permeability in accordance with a Gaussian or 
related stochastic model, coherent with the a priori geologic data, and carrying out, by 
means of a suitable simulator, a simulation of the fluid flows, 

10 b) identifying zones inside the reservoir, calculating the effective permeabilities of 

these zones and, from the simulator results, estimating the corrections to be brought to 
these effective permeabilities to improve calibration in relation to the data, and 

- a second stage comprising : c) propagating the corrections to the whole of the 
grid cells of the permeability field, by means of an iterative optimization process 

15 comprising minimizing a fimctional which quantifies the difference between the 
effective permeabilities required to obtain said calibration and the effective 
permeabilities calculated for the permeability field considered, using a technique of 
gradual deformation of realizations of the stochastic model. 

The zones are defmed either manually or automatically fi-om the flow simulator. 

20 According to an embodiment, simulation of the flows is carried out by means of a 
streamline simulator and the zones of the medium are identified by a set of grid cells 
traversed by one or more streamlines of fixed geometry. 



According to another embodiment, said zones are identified as volume portions on 
the periphery of wells running through said medium, within the framework of well tests. 

According to another embodiment, at least one gradual deformation parameter is 
assigned to each zone. 

5 In relation to the prior art, the method according to the invention allows to calibrate 

a realization of a stochastic model to dynamic data while keeping the coherence in 
relation to the stochastic model and by reducing significantly the number of direct flow 
simulations to be carried out Unlike the approaches developed to date, a deformation of 
the realization does not systematically involve a new flow simulation. It thus affords the 

10 possibility of better exploring the space of the realizations and of rapidly determining 
not only one constrained realization, but several ones. Furthermore, the method 
according to the invention allows to deform the realization from a large number of 
parameters, which accelerates the intermediate optimization. It also allows deformation 
of the realization by zones, and the latter can be the same as the zones predetermined for 

1 5 calculation of the effective permeabilities during the first stage. 

BRIEF DESCRIPTION OF THE FIGURES 

Other features and advantages of the method according to the invention will be 
clear from reading the description hereafter, given by way of non limitative application 
example, with reference to the accompanying drawings wherein : 

20 - Figure 1 shows an example of an optimization scheme in two stages according to the 
method, and 



- Figure 2 shows a comparison between a real fractional flow and the corresponding 
simulated fractional flow. 



DETAILED DESCRIPTION 

Traditionally, calibration of reservoir models is an iterative process wherein a flow 
5 simulation is carried out each time the reservoir is disturbed. For simplicity's sake, we 
consider that the numerical reservoir model comes down to a realization of a stochastic 
model for the permeability. In other words, a flow simulation is required for any 
variation of the permeability fleld. On the contrary, the method according to the 
invention can propose several permeability field variations using a single flow 

10 simulation. It therefore consists of two stages (Fig.l). During the first stage, zones of 
the reservoir are determined and the correction to be brought to the effective 
permeabilities of these zones in order to improve data calibration is estimated. Then, an 
optimization process is started to propagate the disturbance determined for the effective 
permeabilities of these zones to the permeability field representing the reservoir. The 

15 first stage requires a flow simulation because it depends on a comparison between real 
data and the corresponding synthetic responses. On the other hand, tke second stage 
requires no additional flow simulation. Finally, the optimization performed in the 
second stage is parameterized according to the gradual deformation method, which 
allows to preserve the coherence of the permeability values distribution (or any other 

20 property considered) in relation to the spatial variability model. 

The general algorithm is simimarized as follows : 

An initial permeability field in accordance with a stochastic model being generated, 
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a) stage 1 - a flow simulation is performed, the effective permeabilities of these 
zones are calculated, zones are identified in the field considered and the variation to be 
applied to the effective permeabilities of these zones in order to reduce the difference 
between the production data measured in the field and the simulated production data is 

5 estimated; 

b) stage 2 - the disturbance required at the level of the effective permeabilities of 
the zones is propagated to the whole of the permeability field by means of an 
optimization process. The fimction to be minimized quantifies the difference between 
the desired effective permeabilities and the effective permeabilities calculated for the 

10 permeability field considered. The gradual deformation technique is used to modify the 
permeability field; 

c) return to b) as long as the calibration is not satisfactory. 

The effective permeabilities of the zones can be calculated using a simulator or they 
can be obtained by means of averaging techniques well-known to the man skilled in the 
15 art. 

The method according to the invention is adaptable to any flow simulator insofar as 
means to diefine effective permeabilities for different zones of the reservoir are 
available, the latter being identified manually by the user or automatically fixsm a given 
criterion. For example, for a well test simulator, rings of increasing radius, centred on 
20 the wells, can be selected to define the zones. The effective permeabilities can in this 
case be related to the apparent permeabilities. 
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In the example developed hereafter, we focus on the case of a streamline simulator. 
The zones considered here are determined by the streamlines themselves. Each one of 
the aforementioned two st^es is described in detail in the two sections hereunder. 

Stage 1 : Modification of the effective permeabilities 

5 The objective of the first part is to describe the geometry of the streamlines and to 

identify the modification to be applied to the effective permeability of these lines so as 
to better calibrate the production data. A line or a set of streamlines define a zone. 

Modelling flows by streamlines involves certain characteristics. The displacement 
of a fluid along a streamline is a one-dimensional problem; the streamlines do not 

10 communicate with one another. When the mobility ratio of the fluids in presence is one 
and when the boundary conditions do not change, the geometry of the streamlines is 
fixed. When the mobility ratio is different from one, there are two alternatives : the 
geometry of the streamlines is fixed and the flow variations during the fluid 
displacement is allowed or the flow is uniformly distributed between the current tubes 

15 and the geometry of the streamlines is periodically updated. We choose the &st 
configuration straightaway. 

To simplify the problem, we consider groups of streamlines rather than the 
streamlines individually : we thus reduce the number of parameters. The effective 
permeability of a group of streamlines is expressed as the harmonic mean weighted by 
20 the fluid volumes of the permeabilities of the grid cells traversed by the streamlines : 
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Eq.l 



K^^ is the effective permeability of the i-th group of streamlines. NgSL is the 
number of streamlines contained in this i-th group, Nk is the number of grid cells 
traversed by the k-th streamline of the i-th group, qy is the flow for the k-th streamline 
5 at the level of the j-th grid cell. Axkj is the flight time for the k-th streamline through the 
j-th grid cell. 

For a given initial realization, the flow simulation allows to determine the geometry 
of the streamlines and to calculate their effective permeabilities. The effective 
permeability variation that would bring a production data calibration improvement 

10 remains to be evaluated. The fractional flows observed and the corresponding simulated 
fractional flows (Fig.2) are therefore compared for the producing wells. The streamlines 
are then arranged in the increasing order of their breakthrough times and the fractional 
flow curves are discretized. The segments Aq thus defined are associated with groups of 
streamlines. In the absence of accordance between the simulated flows and the real 

15 flows, the effective permeabilities of the groups of streamlines are considered to be 
responsible for the differences. A correction applied to all of the streamline groups 
allows this difference to be reduced. Consider the group of streamlines associated with a 
flow increment Aqi (Fig.2). The effective permeability desired for streamline group i in 
order to improve calibration is : 



20 



i, simulated 



.reference 
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where ti,reference and ti^simuiated are the mean breakthrough times associated with increment 
Aqi for the real and simulated fractional flow curves. The same procedure is carried out 
for each group of streamlines. 

Stage 2 : Intermediate optimization - Modification of the permeability field 

5 The first stage determines a variation to be applied to the effective permeabilities of 

' the streamlines to calibrate the fractional flows. The goal of the second stage is to 
transfer this variation of the streamlines to the permeability field while preserving the 
coherence of the permeability field in relation to the stochastic model. An optimization 
problem is therefore defined. The fiinction referred to as intermediate function is to be 

10 minimized here : 

Ng is the number of streamline groups. K^^^ is the desired effective permeability 
for the group i of streamlines in order to reduce the differences between the measured 
and simulated fractional flows : this value has been determined during the first stage. 
Kf^i„uiated the effective permeability of the group i of streamlines for the permeability 
field considered. This optimization problem is non-linear and can involve as many 
parameters as the permeability field contains grid cells. By integrating the deformation 
method as the parameterization technique, we can reduce the number of parameters and 
provide a permeability field modification that matches the stochastic model defining the 
spatial distribution of the heterogeneities in the permeability field. Furthermore, in order 
to be able to modify the permeability field by zones, the zones corresponding to the 
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various streamline groups, the gradual deformation is not applied to the permeability 
field itself, but to the underlying Gaussian white noise (Fig. 1). 

The simplest version of the gradual deformation method consists in combining two 
Gaussian white noises, z\ and zj, according to the relation : 

z(xXp) = Zj(x)cos(;r/7) + Z2(x)sin(;r>c?) 

where p is the deformation parameter and x the position. The gradient of z with respect 
to p is : 

— = -Tiz, sin(^p) + 71Z2 cos(;r/?). 
dp 

The Gaussian white noise z is then converted to a Gaussian realization y of 
covariance model C, average zero and variance one, by means of a convolution product 
developed at the heart of the FFTMA generator : 

y = f*z. 

f results from the decomposition of the covariance function. The derivative of y 
with respect to the deformation parameter is : 

dp dp' 

Any other geostatistical generator producing Gaussian or related realizations can be 
used in place of the FFTMA generator, provided that calculation of the gradients is 
integrated therein. The FFTMA generator however affords the advantage of fast 
execution, even for realizations discretized on a very large number of grid cells. 



Realization y, which is for the moment centred and reduced, i.e. of average zero and 
variance one, can then be converted to a Gaussian realization w of average m and 
variance ; 

w(x) = m + qv(x). 

5 The derivative with respect to the deformation parameter becomes : 

—{x) = a—{x). 
op dp 

At this stage, static data observed at precise points, in wells for examples, also have 
to be taken into accovmt. This information is in general integrated in the realization 
generated by means of a kriging technique. The constrained realization Wc is deduced 
10 from : 

^c(x) = (x) + (wix) - Wf, (x)) 

where WdK and wk are the realizations obtained, for the first one, from kriging of the 
real observations and, for the second, from kriging of the values of w at the observation 
points. The kriging estimator, in the dual frame, is expressed as follows : 

15 Wf,(x) = Y.p,C{x-Xi) + m. 

C is the covariance function. The Xi are the positions of the n observations. The 
weights Pi are independent of the position, but they depend on the deformation 
parameters. We can show that the derivative of w with respect to the deformation 
parameter is obtained from : 



dp tiSp 



The permeability field k is deduced j&om the lognonnal transformation of Wc : 
k{x) = exp(w^(x)). 

TTie permeability gradient with respect to the deformation parameter is expressed as 
follows : 



op op 



Equation 1) allows to calculate the effective permeability of a group of streamlines 
for a given permeability field. The derivative of the effective permeability with respect 
to the deformation parameter is deduced therefrom : 



dKf ^^^ hh dp 
dp ' ^|, g,.,Ar,,, 



These various relations show how to deform a permeability field and to calculate 
the gradients of the effective permeabilities of the streamlines with respect to a 
deformation parameter. All these relations can be readily generalized in cases where 
several deformation parameters are involved. One can notably decide to assign a 
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deformation parameter to each zone defined for calculation of the effective 
permeabilities. 

Several techniques can be considered for minimizing the intermediate function, but 
since the gradients of the effective permeabilities are available for an insignificant 
5 calculation cost, a Gauss-Newton type approach may be judicious. To determine the 
disturbance to be applied to the deformation parameters so as to reduce the intermediate 
fimction, the following system is solved : 

HLp = -V(F/). 

Ap is the disturbance to be defined, V(FI) comprises the gradients of the 

10 intermediate function with respect to the deformation parameters and H is an 
approached matrix of the Hessian matrix : 

H = G'WG. 

G is the sensitivity matrix : it includes the derivatives of the effective permeabilities 
of the streamline groups with respect to the deformation parameters. W is the weight 
15 matrix : it is here equal to the identity matrix. 

Finally, the algorithm developed to minimize the intermediate function is described 
as follows. Stage 1 allovra to define, on the one hand, the desired effective 
permeabilities for the streamline groups and, on the other hand, zones assigned to these 
groups. 

20 a) At least one deformation parameter is assigned by zone. 
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b) An initial Gaussian white noise and at least one complementary Gaussian white 
noise are defined. 

c) A gradual deformation is carried out by combining the previous Gaussian white 
noises. The derivatives of the Gaussian white noise z resulting from the gradual 

5 deformation with respect to the deformation parameters are simultaneously calculated. 

d) The Gaussian white noise z is converted to a permeability field k and the 
derivatives of the permeabilities with respect to the deformation parameters are also 
determined. 

e) The effective permeabilities of the streamline groups defined in stage 1 and their 
10 derivatives are calculated. 

f) The disturbance to be applied to ihe deformation parameters to reduce the 
intermediate function is then estimated. 

g) At this stage, several alternatives emerge. If the intermediate function is not 
weak enough and if it does not appear to have converged, the deformation parameters 

15 are updated and one returns to c). If the intermediate function is not weak enough but 
seems to have converged, one returns to b), i.e. the initial Gaussian white noise is 
updated and a new complementary Gaussian white noise is randomly selected. Finally, 
if the objective function is sufficiently weak or if it appears that it is not going to 
decrease further, stage 2 comes to an end. 

20 The procedure described here focuses on the case of streamline flow simixlations. 

Streamlines actually appear to be a very natural tool for defining zones. The heart of the 
present invention involving identification of zones and calculation of the effective 
permeabilities for these zones, this choice seems to be logical. Other types of 
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application can however be considered. For example, the information relative to the 
various zones could be provided by well tests. The effective permeability can in this 
case be calculated as a function of the radius of investigation aroimd the well ; the zones 
considered are then rings. The flow simulator used for well test simulation can be a 

5 standard flow simulator. The procedure to be followed is similar to the procedure 
described for the streamlines. Flow simulation allows to identify zones and to determine 
the effective permeability for these zones, which can be compared with the data 
measured in the field. Then, minimization of an intermediate objective function 
according to the approach described above allows to propagate the correction to be 

10 applied to the effective permeabilities to the absolute permeabilities of the grid cells in 
the zones while respecting the a priori spatial variability model. 



